Libraries

We need a bunch of these…

JSON processing

First jsonlite for processing the colours data pulled from the website

library(jsonlite)

Data wrangling

Next, a bunch of data munging packages from the ‘tidyverse’

library(tidyr)
library(plyr) # rbind.fill useful for ragged data with missing entries
library(dplyr)  
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(magrittr)
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:tidyr':
## 
##     extract
library(tibble)
library(stringr)

Spatial

Next, the basic R spatial packages

library(sf)
## Linking to GEOS 3.8.0, GDAL 3.3.1, PROJ 8.1.0
library(raster)
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following object is masked from 'package:magrittr':
## 
##     extract
## The following object is masked from 'package:dplyr':
## 
##     select
## The following object is masked from 'package:tidyr':
## 
##     extract
library(tmap)
library(tmaptools) # for geocode_OSM
tmap_mode("view") # for web maps
## tmap mode set to interactive viewing

Interpolation

Finally, a range of interpolation tools

library(spatstat) # for IDW
## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## spatstat.geom 2.2-2
## 
## Attaching package: 'spatstat.geom'
## The following objects are masked from 'package:raster':
## 
##     area, rotate, shift
## Loading required package: spatstat.core
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:raster':
## 
##     getData
## The following object is masked from 'package:dplyr':
## 
##     collapse
## Loading required package: rpart
## spatstat.core 2.3-0
## Loading required package: spatstat.linnet
## spatstat.linnet 2.3-0
## 
## spatstat 2.2-0       (nickname: 'That's not important right now') 
## For an introduction to spatstat, type 'beginner'
library(maptools)
## Checking rgeos availability: TRUE
library(akima)    # for triangulation
library(fields)   # for splines
## Loading required package: spam
## Loading required package: dotCall64
## Loading required package: grid
## Spam version 2.7-0 (2021-06-25) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
## 
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
## 
##     backsolve, forwardsolve
## Loading required package: viridis
## Loading required package: viridisLite
## See https://github.com/NCAR/Fields for
##  an extensive vignette, other supplements and source code

Getting the colours

See the Dulux website for what this is all about

First I had a poke around on the website to figure out where the colour details were to be found

colour_groups <- c("blues", "browns", "greens", "greys", "oranges",
                  "purples", "reds", "whites-neutrals", "yellows")
base_url <- "https://www.dulux.co.nz/content/duluxnz/home/colour/all-colours.categorycolour.json/all-colours/"

Retrieve and save the colours

The loop on the next slide

  • steps through the group names
  • retrieves the relevant JSON file
  • writes it out locally
  • adds the colours information to a list
  • then we make the list into a single table using bind_rows
colours <- list()
for (i in 1:length(colour_groups)) {
  colour_group <- colour_groups[i]
  json_url <- str_c(base_url, colour_group)
  json <- fromJSON(json_url, flatten = TRUE)
  # make a local copy (just for convenience)
  json_file_name <- str_c(colour_group, ".json")
  write_json(json, json_file_name)
  # get the colours information and add to list
  the_colours <- rbind.fill(json$categoryColours$masterColour.colours)
  colours[[i]] <- the_colours
  Sys.sleep(0.5) # pause to not annoy the the server  
}
df_colours <- bind_rows(colours)
write.csv(df_colours, "dulux-colours-raw.csv", row.names = FALSE)

Check we’re all good

head(df_colours)
##       id red green blue lrv      baseId           name woodType coats
## 1 149253 205   210  206  67 vivid_white Pukaki Quarter     None    NA
## 2 149254 220   240  242  86 vivid_white      Canoe Bay     None    NA
## 3 149255 226   240  245  87 vivid_white      Mt Dobson     None    NA
## 4 149256 220   230  235  80 vivid_white        Raetihi     None    NA
## 5 149257 217   219  223  74 vivid_white   Taiaroa Head     None    NA
## 6 149258 180   200  219  60 vivid_white   Gulf Harbour     None    NA

Tidying up the names

There are paint names with modifiers as suffixes for different shades of particular colours, and we need to handle this

The modifiers are

paint_modifiers <- c("Half", "Quarter", "Double")

A tidyverse pipeline

Here’s one way to clean this up (there are others…)

df_colours_tidied <- df_colours %>%
  ## separate the name components, filling from the left with NAs if <5
  separate(name, into = c("p1", "p2", "p3", "p4", "p5"), sep = " ",
           remove = FALSE, fill = "left") %>%
  ## replace any NAs with an empty string
  mutate(p1 = str_replace_na(p1, ""),
         p2 = str_replace_na(p2, ""),
         p3 = str_replace_na(p3, ""),
         p4 = str_replace_na(p4, "")) %>%
  ## if p5 is a paint modifiers, then recompose name
  ## from p1:p4 else from p1:p5
  ## similarly keep modifier where it exists
  mutate(placename = if_else(p5 %in% paint_modifiers,
                       str_trim(str_c(p1, p2, p3, p4, sep = " ")),
                       str_trim(str_c(p1, p2, p3, p4, p5, sep = " "))),
         modifier = if_else(p5 %in% paint_modifiers,
                       p5, "")) %>%
  ## remove some places that are tricky to deal with later
  filter(!placename %in% c("Chatham Islands",
                           "Passage Rock",
                           "Auckland Islands",
                           "Cossack Rock")) %>%
  ## throw away variables we no longer and reorder
  select(name, placename, modifier, red, green, blue)

# save it so we have it for later
write.csv(df_colours_tidied, "dulux-colours.csv", row.names = FALSE)

Build the spatial dataset

Add x and y columns to our data for the coordinates

df_colours_tidied_xy <- df_colours_tidied %>%
  mutate(x = 0, y = 0)

Geocode with tmaptools::geocode_OSM

Code on the next slide

  • goes through all the unique placenames
  • appends as many x y coordinates as we have space for (due to the modifiers) from the geocoding results

Best not to re-run this (it takes a good 10 minutes and it’s not good to repeatedly geocode and hit the OSM server)

for (placename in unique(df_colours_tidied_xy$placename)) {
  address <- str_c(placename, "New Zealand", sep = ", ")
  geocode <- geocode_OSM(address, as.data.frame = TRUE, return.first.only = FALSE)
  num_geocodes <- nrow(geocode)
  matching_rows <- which(df_colours_tidied_xy$placename == placename)
  for (i in 1:length(matching_rows)) {
    if (!is.null(geocode)) {
      if (num_geocodes >= i) {
        df_colours_tidied_xy[matching_rows[i], ]$x <- geocode$lon[i]
        df_colours_tidied_xy[matching_rows[i], ]$y <- geocode$lat[i]
      }  
    }
  }
  Sys.sleep(0.5) # so as not to over-tax the geocoder
}
# Remove any we missed
df_colours_tidied_xy <- df_colours_tidied_xy %>%
  filter(x != 0 & y != 0)

write.csv(df_colours_tidied_xy, "dulux-colours-xy.csv", row.names = FALSE)

Making maps

Make the dataframe into a sf point dataset

dulux_colours_sf <- st_as_sf(df_colours_tidied_xy,
                     coords = c("x", "y"), # columns with the coordinates
                     crs = 4326) %>%       # EPSG:4326 for lng-lat
  st_transform(2193) %>% # convert to NZTM
  ## and make an RGB column
  mutate(rgb = rgb(red / 255, green / 255, blue/ 255))

# jitter any duplicate locations
duplicate_pts <- which(duplicated(dulux_colours_sf$geometry) |
                       duplicated(dulux_colours_sf$geometry, fromLast = TRUE))
jittered_pts <- dulux_colours_sf %>%
  slice(duplicate_pts) %>%
  st_jitter(50)
dulux_colours_sf[duplicate_pts, ]$geometry <- jittered_pts$geometry

st_write(dulux_colours_sf, "dulux-colours-pts.gpkg", delete_dsn = TRUE)

Remember to update the dataframe of points with the jittered points

jittered_pts <- dulux_colours_sf %>%
  st_coordinates() %>%
  as_tibble()
df_colours_tidied_xy <- df_colours_tidied_xy %>%
  mutate(x = jittered_pts$X, y = jittered_pts$Y)
write.csv(df_colours_tidied_xy, "dulux-colours-xy-jit.csv", row.names = TRUE)

And at last a map!

nz <- st_read("nz.gpkg")
## Reading layer `nz' from data source 
##   `/home/osullid3/Documents/code/dulux-colours-map/nz.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 0 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 1089972 ymin: 4748123 xmax: 2091863 ymax: 6223160
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
## Reading layer `dulux-colours-pts' from data source 
##   `/home/osullid3/Documents/code/dulux-colours-map/dulux-colours-pts.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 904 features and 7 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1094875 ymin: 4759975 xmax: 2073657 ymax: 6191716
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
tm_shape(nz) +
  tm_borders() +
  tm_shape(dulux_colours_sf) +
  tm_dots(col = "rgb")

We can do better (depending on taste!)

Points aren’t really much fun, instead:

Voronoi polygons

We can make Voronois and clip to NZ

dulux_colours_vor <- dulux_colours_sf %>%
  st_union() %>%
  st_voronoi() %>%
  st_cast() %>%
  st_as_sf() %>%
  st_join(dulux_colours_sf, left = FALSE) %>%
  st_intersection(st_read("nz.gpkg"))

st_write(dulux_colours_vor, "dulux-colours-vor.gpkg", delete_dsn = TRUE)

And map

## Reading layer `dulux-colours-vor' from data source 
##   `/home/osullid3/Documents/code/dulux-colours-map/dulux-colours-vor.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 902 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 1089972 ymin: 4748123 xmax: 2091863 ymax: 6223160
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
tm_shape(dulux_colours_vor) +
  tm_polygons(alpha = 0, border.col = "lightgrey", lwd = 0.5) +
  tm_shape(dulux_colours_vor) +
  tm_polygons(col = "rgb", id = "placename", 
              alpha = 0.75, border.col = "grey", lwd = 0.2)

Triangulation

Using the akima::interp package

components = c("red", "green", "blue")

layers = list()
for (component in components) {
  # the dimensions, nx, ny give ~500m resolution
  layers[[component]] <- raster(
    interp(df_colours_tidied_xy$x, df_colours_tidied_xy$y, 
           df_colours_tidied_xy[[component]],
           nx = 2010, ny = 2955, linear = TRUE))
}
rgb.t <- brick(layers)
crs(rgb.t) <- st_crs(nz)$wkt
rgb.t <- mask(rgb.t, nz)

writeRaster(rgb.t, "dulux-colours-tri.tif", overwrite = TRUE)

And map

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition
tm_shape(dulux_colours_vor) +
  tm_polygons(alpha = 0, border.col = "lightgrey", lwd = 0.5) +
  tm_shape(rgb.t) + 
  tm_rgb()
## stars object downsampled to 825 by 1213 cells. See tm_shape manual (argument raster.downsample)

IDW

Here we use the spatstat::idw function

We have to make spatstat::ppp (point pattern) objects

## We need a window for the ppps
W <- nz %>%
  st_geometry() %>%
  st_buffer(1000) %>%
  st_union() %>%
  as("Spatial") %>%
  as.owin()
 
layers <- list()
for (component in components) {
  pp <- ppp(x = df_colours_tidied_xy$x, 
            y = df_colours_tidied_xy$y, 
            window = W, 
            marks = df_colours_tidied_xy[[component]])
  ## eps is the approximate resolution
  layers[[component]] <- raster(idw(pp, eps = 2500, power = 4))
}
rgb.idw <- brick(layers)
crs(rgb.idw) <- st_crs(nz)$wkt

writeRaster(rgb.idw, "dulux-colours-idw.tif", overwrite = TRUE)

And map

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition
tm_shape(dulux_colours_vor) +
  tm_polygons(alpha = 0, border.col = "lightgrey", lwd = 0.5) +
  tm_shape(rgb.idw) + 
  tm_rgb()

Splines

For this we use the fields::Tps (thin plate spline) function

We need an empty raster to interpolate into

r <- nz %>% 
  st_make_grid(cellsize = 2500, what = "centers") %>%
  st_coordinates() %>%
  as_tibble() %>%
  rename(x = X, y = Y) %>%
  mutate(z = 0) %>%
  rasterFromXYZ()

Then we can interpolate as before

layers <- list()
for (component in components) {
  spline <- Tps(df_colours_tidied_xy[, c("x", "y")], 
                df_colours_tidied_xy[[component]], 
                scale.type = "unscaled", m = 3)
  layers[[component]] <- interpolate(r, spline)
}
rgb.s <- brick(layers)
crs(rgb.s) <- st_crs(nz)$wkt
rgb.s <- mask(rgb.s, nz)

writeRaster(rgb.s, "dulux-colours-spline.tif", overwrite = TRUE)

And map

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition
tm_shape(dulux_colours_vor) +
  tm_polygons(alpha = 0, border.col = "lightgrey", lwd = 0.5) +
  tm_shape(rgb.s) + 
  tm_rgb()
## Warning: Raster values found that are outside the range [0, 255]

Areal interpolation

We can also do weighted mixing by areas of overlap with any other set of areas from the Voronois.

sa2 <- st_read("sa2-generalised.gpkg")

dulux_colours_sa2 <- dulux_colours_vor %>%
  select(red:blue) %>%
  st_interpolate_aw(sa2, extensive = FALSE) %>%
  mutate(rgb = rgb(red / 256, green / 256, blue / 256),
         name = sa2$SA22018_V1_00_NAME)

st_write(dulux_colours_sa2, "dulux-colours-sa2.gpkg", delete_dsn = TRUE)

And map

## Reading layer `dulux-colours-sa2' from data source 
##   `/home/osullid3/Documents/code/dulux-colours-map/dulux-colours-sa2.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 2154 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 1089972 ymin: 4748123 xmax: 2091863 ymax: 6223160
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
tm_shape(dulux_colours_vor) +
  tm_fill(col = "rgb", id = "placename") + 
  tm_shape(dulux_colours_sa2) +
  tm_polygons(col = "rgb", id = "name", border.col = "white", lwd = 0.2)

Credits and more